Comparison between MAP calibrated model outcomes and calibration targets - not time dependent

Published

18/04/2023 T12:02

Code
library(here)

source(here("02_scripts/03b_fun_calib_ntd.R"))
theme_set(theme_minimal())

opt_params_notimedepend <- readRDS(here("01_data/map_calib_ntd.RDS"))

Used overall deaths, prevalence of rx opioid use, and total OAT as calibration targets

Code
mod_opt_map_ntd <- mod_calib_ntd(as.numeric(opt_params_notimedepend$par))

1 Prevalence of Rx opioids use

Code
# prevalence of people on prescribed opioids

prop_opioids_rx_calib <- prop_opioids_rx_uncalib <- data.frame(matrix(NA, byrow = T,
                                             nrow = 12,
                                             ncol = length(2015:2018),
                                             list(NULL,
                                                  paste0("prop_opioids_rx_",
                                                         str_sub(2015:2018, 3, 4)))))
tot_pop_calib_tbl <- tot_pop_uncalib_tbl <- data.frame(matrix(NA, byrow = T,
                                         nrow = 12,
                                         ncol = length(2015:2018),
                                 list(NULL,
                                      paste0("tot_pop_",
                                             str_sub(2015:2018, 3, 4)))))

for (i in 1:length(2015:2018)) {
  yr <- 2014 + i
  prop_opioids_rx_uncalib[, i] <- assign(
    paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
    rowSums(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
                             c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
                               "BPO_OTHER", "BPO_PALLIATIVE")]) /
      rowSums(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
    
  )
  
  tot_pop_uncalib_tbl[, i] <- assign(
    paste0("tot_pop_", str_sub(yr, 3, 4)),
    rowSums(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
  )
  
  prop_opioids_rx_calib[, i] <- assign(
    paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
    rowSums(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
                             c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
                               "BPO_OTHER", "BPO_PALLIATIVE")]) /
      rowSums(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
    
  )
  
  tot_pop_calib_tbl[, i] <- assign(
    paste0("tot_pop_", str_sub(yr, 3, 4)),
    rowSums(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
  )
  }

prop_opioids_rx_uncalib_tbl <- prop_opioids_rx_uncalib %>% 
  pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>% 
  arrange(grp) %>% 
  mutate(year = rep(2015:2018,
                              each = 12),
         year_mon = paste(year, month.abb[1:12], sep = "_")) %>% 
  bind_cols(., tot_pop_uncalib_tbl %>% 
              pivot_longer(1:4, names_to = "tot_pop",
                           values_to = "tot_pop_val") %>% 
              arrange(tot_pop) %>% select(-tot_pop))


prop_opioids_rx_calib_tbl <- prop_opioids_rx_calib %>% 
  pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>% 
  arrange(grp) %>% 
  mutate(year = rep(2015:2018,
                              each = 12),
         year_mon = paste(year, month.abb[1:12], sep = "_")) %>% 
  bind_cols(., tot_pop_calib_tbl %>% 
              pivot_longer(1:4, names_to = "tot_pop",
                           values_to = "tot_pop_val") %>% 
              arrange(tot_pop) %>% select(-tot_pop))


prop_opioids_rx_target_tbl <- calib_target_tbl %>% 
              filter(group == "prev_on_opioidsrx") %>% 
              select(year, target, group) %>% 
              rename(grp = group,
                     target_val = target) %>% 
              mutate(year_mon = paste(year, month.abb[6], sep = "_"))

prop_opioids_rx_uncalib_wei_mean <- prop_opioids_rx_uncalib_tbl %>% 
  group_by(year) %>% 
  summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>% 
  ungroup() %>% 
  mutate(grp = "Uncalibrated Model") %>% 
  bind_rows(., prop_opioids_rx_target_tbl %>% 
              select(year, target_val) %>% 
              mutate(grp = "target")) %>% 
  bind_rows(., prop_opioids_rx_calib_tbl %>% 
              group_by(year) %>% 
              summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>% 
              ungroup() %>% 
              mutate(grp = "Calibrated Model"))

prop_opioids_rx_target_tbl$year_mon <- factor(prop_opioids_rx_target_tbl$year_mon,
                                              levels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"),
                                              labels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"))
prop_opioids_rx_uncalib_tbl$year_mon <- factor(prop_opioids_rx_uncalib_tbl$year_mon,
                                              levels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"),
                                              labels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"))
prop_opioids_rx_calib_tbl$year_mon <- factor(prop_opioids_rx_calib_tbl$year_mon,
                                              levels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"),
                                              labels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"))

prop_opioids_rx_tbl <- prop_opioids_rx_uncalib_tbl %>% 
  mutate(group = "Uncalibrated Model") %>% 
  bind_rows(.,prop_opioids_rx_calib_tbl %>% 
              mutate(group = "Calibrated Model"))

wprop_opioids_rx_tbl <- prop_opioids_rx_target_tbl %>% 
  mutate(grp = "Target") %>% 
  bind_rows(., prop_opioids_rx_uncalib_tbl %>% 
               group_by(year) %>% 
               summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>% 
               ungroup() %>% 
               mutate(year_mon = prop_opioids_rx_target_tbl$year_mon,
                      grp = "Uncalibrated Model")) %>% 
  bind_rows(., prop_opioids_rx_calib_tbl %>% 
               group_by(year) %>% 
               summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>% 
               ungroup() %>% 
               mutate(year_mon = prop_opioids_rx_target_tbl$year_mon,
                      grp = "Calibrated Model"))

prop_opioids_rx_tbl$group <- factor(prop_opioids_rx_tbl$group,
                                  labels = c("Calibrated Model", "Uncalibrated Model"),
                                  levels = c("Calibrated Model", "Uncalibrated Model"))


p1 <- ggplot() +
  geom_point(data = prop_opioids_rx_tbl, aes(x = year_mon, y = target_val,
                                             color = group, fill = group)) +
  geom_point(data = wprop_opioids_rx_tbl %>% filter(grp == "Target"),
             aes(x = year_mon, y = target_val, color = grp, fill = grp), size = 2) +
   scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  xlab("Year_Month") + ylab("Prevalence of Rx opioid use") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

plotly::ggplotly(p1)

1.0.1 Comparison between uncalibrated annual weighted average and observed targets of prevalence of Rx opioids use

Code
wprop_opioids_rx_tbl$grp <- factor(wprop_opioids_rx_tbl$grp,
                                   labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
                                   levels = c("Target", "Calibrated Model", "Uncalibrated Model"))

p2 <- ggplot() +
 geom_point(data = wprop_opioids_rx_tbl,
             aes(x = year_mon, y = target_val,
                 color = grp, fill = factor(grp)), size = 2) +
  xlab("Year") + ylab("Prevalence of Rx opioid use") +
  scale_color_brewer(palette = "Dark2", name = "") +
  scale_fill_brewer(palette = "Dark2", name = "")
  # ylim(0,0.2) +

plotly::ggplotly(p2)

2 Deaths and OD-Deaths

Code
################ Deaths 
# number of deaths

num_deaths_calib <- num_deaths_uncalib <- rep(NA, length(2016:2020))
for (i in 1:length(2016:2020)){
  yr <- 2015 + i
  num_deaths_uncalib[i] <- mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                       year_mon_cycle_tbl$mon == 12] + 1,
                                            "BO_DEATH"] -
    mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_DEATH"]
  
  num_deaths_calib[i] <- mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                       year_mon_cycle_tbl$mon == 12] + 1,
                                            "BO_DEATH"] -
    mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_DEATH"]
  }
  
# number of OD-deaths
  
num_od_deaths_calib <- num_od_deaths_uncalib <- rep(NA, length(2016:2021))
for (i in 1:length(2016:2021)){
  yr <- 2015 + i 
  num_od_deaths_uncalib[i] <- mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
  num_od_deaths_calib[i] <- mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
}

deaths_target_uncalib_tbl <- calib_target_tbl %>% 
  select(year, target, group) %>% 
  rename(target_val = target,
         target = group) %>% 
  filter(target %in% c("total_deaths", "total_od_deaths")) %>% 
  mutate(target = ifelse(target == "total_deaths", "Total deaths",
                         ifelse(target == "total_od_deaths",
                                "Total OD-deaths", NA)),
         group = "Target") %>% 
  bind_rows(., data.frame(target_val = num_deaths_uncalib) %>% 
              mutate(target = "Total deaths") %>% 
              bind_rows(., data.frame(target_val = num_od_deaths_uncalib) %>%
                          mutate(target = "Total OD-deaths")) %>% 
              mutate(year = c(2016:2020, 2016:2021),
                     group = "Uncalibrated Model")) %>% 
   bind_rows(., data.frame(target_val = num_deaths_calib) %>% 
              mutate(target = "Total deaths") %>% 
              bind_rows(., data.frame(target_val = num_od_deaths_calib) %>%
                          mutate(target = "Total OD-deaths")) %>% 
              mutate(year = c(2016:2020, 2016:2021),
                     group = "Calibrated Model"))
  
deaths_target_uncalib_tbl$group <- factor(deaths_target_uncalib_tbl$group,
                                   labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
                                   levels = c("Target", "Calibrated Model", "Uncalibrated Model"))

p3 <- ggplot() +
  geom_point(data = deaths_target_uncalib_tbl %>% filter(target == "Total deaths"),
             aes(x = year, y = target_val, color = factor(group),
                 fill = factor(group))) +
  xlab("Year") + ylab("Total Deaths") +
  scale_color_brewer(palette = "Dark2", name = "") +
  scale_fill_brewer(palette = "Dark2", name = "")
# +
#   facet_wrap(~target, scales = "free") + theme(legend.position="bottom")



plotly::ggplotly(p3)
Code
p4 <- ggplot() +
  geom_point(data = deaths_target_uncalib_tbl %>% filter(target == "Total OD-deaths"),
             aes(x = year, y = target_val, color = factor(group),
                 fill = factor(group))) +
  xlab("Year") + ylab("Total Opioid-related Overdose Deaths") +
  scale_color_brewer(palette = "Dark2", name = "") +
  scale_fill_brewer(palette = "Dark2", name = "")


plotly::ggplotly(p4)

3 OAT

Code
# number of oat
  
num_oat_calib <- num_oat_uncalib <- rep(NA, length(2018:2021))
for (i in 1:length(2018:2021)){
  yr <- 2017 + i 
    num_oat_uncalib[i] <- sum(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                        year_mon_cycle_tbl$mon == 6] + 1,
                             c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
    
    num_oat_calib[i] <- sum(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                        year_mon_cycle_tbl$mon == 6] + 1,
                             c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
}


oat_target_uncalib_tbl <- calib_target_tbl %>% 
  select(year, target, group) %>% 
  rename(target_val = target,
         target = group) %>% 
  filter(target %in% c("total_oat")) %>% 
  mutate(target = ifelse(target == "total_oat",
                                "Total OAT", NA),
         group = "Target") %>% 
  bind_rows(., data.frame(target_val = num_oat_uncalib) %>%
                          mutate(target = "Total OAT",
                                 year = c(2018:2021),
                                 group = "Uncalibrated Model")) %>% 
   bind_rows(., data.frame(target_val = num_oat_calib) %>%
                          mutate(target = "Total OAT",
                                 year = c(2018:2021),
                                 group = "Calibrated Model"))

oat_target_uncalib_tbl$group <- factor(oat_target_uncalib_tbl$group,
                                   labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
                                   levels = c("Target", "Calibrated Model", "Uncalibrated Model"))

p5 <- ggplot() +
  geom_point(data = oat_target_uncalib_tbl,
             aes(x = year, y = target_val, color = factor(group),
                 fill = factor(group))) +
  xlab("Year") + ylab("Total OAT counts") +
  scale_color_brewer(palette = "Dark2", name = "") +
  scale_fill_brewer(palette = "Dark2", name = "")
# +
#   facet_wrap(~target, scales = "free") + theme(legend.position="bottom")



plotly::ggplotly(p5)
Code
cbind.data.frame(parameter = calib_param_tbl_ntd$var_params,
                 basevalue = calib_param_tbl_ntd$basevalue,
                 opt_map = opt_params_notimedepend$par) %>% 
  filter(basevalue != opt_map)
                      parameter  basevalue    opt_map
1          p__BN_PN__BPO_MISUSE 0.00032511 0.00220011
2      p__BPO_CHRONIC__BO_OD_RX 0.00097898 0.00160398
3   p__BPO_PALLIATIVE__BO_DEATH 0.37500000 0.38250000
4 p__BS_OAT_MAINT__BS_OAT_MAINT 0.70000000 0.77000000